#-------------------------------------------------------------
# Script Replicating Results in Extremists Not on Board
# BCK Egerod and H Tran
#
#-------------------------------------------------------------

# Replicating results in the main paper

# load relevant libraries
library(haven); library(readr); library(readxl)
library(margins); library(broom); library(dplyr)
library(stargazer); library(ggplot2); library(cowplot)
library(interflex);library(gridExtra); library(tidyr)
library(specr); library(purrr); library(lfe)
library(here)

# uncomment and run this if you don't have zamfluence installed
 # devtools::install_github("https://github.com/rgiordan/zaminfluence/",
 #                          ref="master",
 #                          subdir="zaminfluence",
 #                          force=TRUE)

library(zaminfluence)

this_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)

sink(paste(this_dir,"EgerodTranLog.txt", sep="/"))

# load data
pol <- read_csv(here(this_dir,"data","BoardServicePolCharacteristics_repli.csv"))

#---------------------
# Figure 1

# obtain AME from the bivariate model -- print in graph
biv_mod <- glm(seats_t_ind ~ 
                 abs_nom, 
               data = pol, family = "binomial")
biv_mar <- summary(margins(biv_mod))

#create graph
ideo<- ggplot(pol, aes(x = dw_nom, y = seats_t_ind)) +
  geom_smooth(se = F, colour = "black") +
  theme_classic() +
  labs(x = "DW-NOMINATE\n(Negative = Liberal; Positive = Conservative)", y = "Probability of Board Service") +
  annotate(geom = "text", x = 0.8, y = 0.15,
           label = paste(round(biv_mar$AME,digits =2),
                         paste("(",round(biv_mar$SE,digits =2), ")", sep =""),
                         sep = "\n"),
           size = 4)

# add marginal distribution of DW-NOMINATE
nom_xhist <- axis_canvas(ideo, axis = "x") +
  geom_histogram(data = pol, aes(x = dw_nom), 
                 fill = "grey")  

ideo <- insert_xaxis_grob(ideo, nom_xhist, position = "bottom")

# draw figure
ggdraw(ideo)

# uncomment to export
# ggsave(plot = ideo,
#        filename = "ProbBoardIdeoDist.eps",
#        device = cairo_ps,
#        width = 6.5, height = 5)


#------------------
# Table 1

# model with immediate seat as DV
mod <- glm(seats_t_ind ~ 
             abs_nom + # extremity
             # expertise
             les + n_bills +
             chair + sub_chair + av_hhi +
             #visibility + connex
             seniority +
             power2 +
             av_connect +
             av_fund +
             #other factors and controls
             state_leg + state_leg_prof + 
             + lost + chamber +  party_txt +  
             majority + female + african_american +
             latinx + 
             vote_share +
             state_po, 
           data = pol, family = "binomial")

#seat within 5 years as DV
mod_5 <- glm(seat_next_5yr ~ 
               abs_nom + # extremity
               les + n_bills +
               chair + sub_chair + av_hhi +
               seniority +
               power2 +
               av_connect +
               av_fund +
               state_leg + state_leg_prof + 
               
               + lost + chamber +  party_txt +  
               majority + female + african_american +
               latinx + 
               vote_share +
               state_po, 
             data = pol, family = "binomial")

# calculate AMEs to put in table alongside logit coefficients
mar_mod <- margins(mod) 
mar_mod <- summary(mar_mod)

order_names <- rownames(as.data.frame(coef(mod)))[-1]

mar_mod_ordered <- mar_mod[match(order_names, mar_mod$factor),]

SEs <- c(sqrt(diag(vcov(mod)))[-1],
         mar_mod_ordered$SE)

# AME for 5-year
mar_mod5 <- margins(mod_5) 
mar_mod5 <- summary(mar_mod5)

order_names5 <- rownames(as.data.frame(coef(mod_5)))[-1]

mar_mod_5_ordered <- mar_mod5[match(order_names5, mar_mod5$factor),]

SEs5 <- c(sqrt(diag(vcov(mod_5)))[-1],
          mar_mod_5_ordered$SE)

# Table 1 containing logit coefficients and AMEs
stargazer(mod, mod,mod_5,mod_5, 
          keep = "nom",
          coef = list(as.vector(coef(mod)),
                      as.vector(c(0,mar_mod_ordered$AME)),
                      as.vector(coef(mod_5)),
                      as.vector(c(0,mar_mod_5_ordered$AME))),
          se = list(as.vector(sqrt(diag(vcov(mod)))),
                    as.vector(c(0,mar_mod_ordered$SE)),
                    as.vector(sqrt(diag(vcov(mod_5)))),
                    as.vector(c(0,mar_mod_5_ordered$SE))),
          omit = c("Constant", "state_po", "factor", "female", 
                   "african_american", "latinx", "party","chamber"),
          covariate.labels = c("Ideological Partisanship",
                               "Legislative Effectiveness Score",
                               "ln Number of Bills", 
                               "Committee Chair", "Subcommittee Chair",
                               "Bill Topic Specialization", "ln Years in Congress",
                               "On Power Committee", "Cosponsor Centrality",
                               "ln Average Campaign Funds",
                               "State Legislature", "State Legislature Professionalism",
                               "Lost Re-Election", "Proportion in Majority", "Average Vote Share"),
          dep.var.labels = c("Immediate Seat", "Seat Within 5 Years"),
          column.labels = c("Log-Odds", "AME", "Log-Odds", "AME"),
          title = "Politician Characteristics and Board Service",
          label = "tab:res",
          type="text", # comment out to get latex table
          #out = "MainTable.tex", # uncomment to export table
          add.lines = list(c("State Fixed Effects?", "Yes", "Yes", "Yes", "Yes"),
                           c("Controls?", "Yes", "Yes", "Yes", "Yes")))


#------------------------
# Figure 2


subs <- data.frame(mod = seq(0.1, 0.3, 0.05),
                   ext = seq(0.7, 0.9, 0.05))

est <- list()

for(i in 1:nrow(subs)){
  
  pol <- pol %>%
    mutate(cat_nom = case_when(abs_nom < quantile(abs_nom, probs = subs[i,1], na.rm = T) ~ 0,
                               abs_nom > quantile(abs_nom, probs = subs[i,2], na.rm = T) ~ 1,
                               TRUE ~ -1))
  
  est_1 <- glm(seats_t_ind ~ factor(cat_nom) +
                 les + n_bills +
                 chair + sub_chair + av_hhi +
                 #visibility + connex
                 seniority +
                 power2 +
                 av_connect +
                 av_fund +
                 #other factors and controls
                 state_leg + state_leg_prof + 
                 
                 + lost + chamber +  party_txt +  
                 majority + female + african_american +
                 latinx + 
                 vote_share +
                 state_po,
               data = pol, family = "binomial")
  
  
  est_1 <- tidy(est_1)[1:3,]
  
  est[[i]] <- rbind(data.frame(est_1, DV = "Immediate", ptile = c(subs[i, 1],subs[i, 1], subs[i, 2]))) 
  
  
}

est <- do.call("rbind", est)
est$spec <- c("C: Loyalists", "A: Moderates", "B: Extremists")


p_ptile <- est %>% filter(spec != "C: Loyalists") %>%
  ggplot(aes(x = ptile, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - 1.68*std.error,
                    ymax = estimate + 1.68*std.error), width = 0) +
  facet_wrap(~ spec, scales = "free_x") +
  geom_hline(yintercept = 0, lty = 3) +
  theme_classic() +
  labs(x = "Percentile Defining Moderates and Extremists\n(Moderates under, Extremists over)",
       y = "Logit Coefficient")
p_ptile

# uncomment to export
# ggsave(p_ptile, 
#        filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/Effect_in_tails.eps",
#        width = 8, height =5,
#        device = cairo_ps)


#-----------------------
# Figure 3

# create identifiers of placebo congresses (pre-HLOGA)
pol$f_hloga <- ifelse(pol$cong>106&pol$cong<=109,1,0)

pol <- pol %>%
  mutate(f_hloga = case_when(cong %in% c(106, 107) ~ 0,
                             cong %in% c(108, 109) ~ 1))

# create placebo treatment
pol$f_did <- pol$sen*pol$f_hloga


# estimate the real (non-placebo) HLOGA DiD
pol <- as.data.frame(pol) #convert to data.frame for interflex
set.seed(888) # set seed

# estimate
ekstreme <- interflex(estimator = "binning",
                      data = pol, 
                      Y = "seats_t_ind",
                      D = "did", X = "abs_nom",
                      Z = c("sen", "hloga"),vartype="bootstrap",
                      cores = 8,nboots=2000,
                      cl="state_po",
                      cutoffs = quantile(pol$abs_nom, probs = c(0.25,0.5,0.75),na.rm=T),#cutoffs = c(0.25,0.5,0.75),
                      na.rm=T)

# estimate placebo (pre-HLOGA) DiD
set.seed(21)
sub_pol <- filter(pol, hloga !=  1) #subset to exclude treated period
placebo <- interflex(estimator = "binning",
                     data = sub_pol, 
                     Y = "seats_t_ind",
                     D = "f_did", X = "abs_nom",
                     Z = c("sen", "f_hloga"),vartype="bootstrap",
                     cores = 8,nboots=2000,
                     cl="state_po",
                     #FE = "cong",
                     cutoffs = quantile(pol$abs_nom, probs = c(0.25,0.5,0.75),na.rm=T),#cutoffs = c(0.25,0.5,0.75),
                     na.rm=T)

# extract estimates for graph
df_ex <- as.data.frame(ekstreme$est.bin$`1`)
df_ex$spec <- "HLOGA"
df_ex$adj_x0 <- df_ex$x0 + c(0.025, 0.03, 0.025, 0)

df_plac <- as.data.frame(placebo$est.bin$`1`)
df_plac$spec <- "Placebo"
df_plac$adj_x0 <- df_plac$x0

df_both <- rbind(df_ex, df_plac)

# make figure 
att_nom <- ggplot(df_both,
                  aes(x = adj_x0, y = coef, colour = spec, group = spec)) +
  geom_point(position = position_dodge(width = 0.1)) +
  geom_errorbar(aes(ymin = CI.lower,
                    ymax = CI.upper), width = 0,
                position = position_dodge(width = 0.1)) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Ideological Partisanship\n(Absolute Distance to Chamber-Period Median of DW-NOMINATE)",
       y = "ATT of HLOGA on Board Service") +
  scale_colour_manual(values = c("black", "grey"),
                      name = " ") +
  scale_x_continuous(limits = c(0,1))

# marginal distribution
nom_xhist <- axis_canvas(att_nom, axis = "x") +
  geom_histogram(data = pol, aes(x = abs_nom), 
                 fill = "grey") 
# 
att_nom <- insert_xaxis_grob(att_nom, nom_xhist, position = "bottom")
ggdraw(att_nom)

# uncomment to export
# ggsave("C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/extreme2.eps",
#         att_nom, device = cairo_ps,
#        height = 4, width = 6)


#-------------------------------------------------------------
# ONLINE APPENDIX for Extremists Not on Board
#

#----------------
# Figure with distribution of DW-NOMINATE scores

p_nom <- ggplot(pol, aes(x = dw_nom)) +
  geom_histogram(fill = "grey") +
  theme_classic() +
  labs(x = "DW-NOMINATE\n(From Liberal to Conservative)",
       y = "Frequency")

p_nom

# uncomment to export
# ggsave(plot = p_nom,
#        filename = "NomDist.eps",
#        device = cairo_ps,
#        width = 6.5, height = 5)


#----------
# Table with descriptive statistics

desc_dat<- select(pol, board_service, seats_t_ind,seat_next_5yr, abs_nom , majority , female , african_american ,
                  latinx , 
                  vote_share , chair , sub_chair , seniority ,
                  state_leg , state_leg_prof , 
                  n_bills , 
                  les , lost , chamber , av_hhi , party_txt , power2 , 
                  av_fund, av_connect)

stargazer(na.omit(desc_dat),
          type = "text", # comment out to create latex table
          covariate.labels = c("Ever on Board?",
                               "Immediate Board Seat?",
                               "Board Seat within 5 Years?",
                               "Ideological Partisanship",
                               "Prop. of Career in Majority",
                               "Female", "African-American", "Latinx",
                               "Average Vote Share", "Committee Chair?",
                               "Sub-Committee Chair?", "Years in Congress (log)",
                               "State Legislature?", "State Legislature Professionalism",
                               "Average Bills Sponsored (log)",
                               "Legislative Effectiveness Scores","Lost Election",
                               "Topic Specialization",
                               "On Power Committee?", "Average Funds Raised",
                               "Cosponsor Centrality"),
          title = "Descriptive Statistics", label = "desc")
          #out = "Descriptives.tex") # uncomment to export table


#----------------------
# Table with correlation matrix


desc_dat<- select(pol, abs_nom , majority , female , african_american ,
                  latinx , 
                  vote_share , chair , sub_chair , seniority ,
                  state_leg , state_leg_prof , 
                  n_bills , 
                  les , lost ,  av_hhi ,  power2 , 
                  av_fund, av_connect)

cor_mat <- cor(desc_dat, use = "pairwise.complete.obs")
cor_mat[upper.tri(cor_mat)] <- NA
rownames(cor_mat) <- c("Ideological Partisanship", 
                       "Prop. of Career in Majority", 
                       "Female", "African-American", "Latinx", 
                       "Average Vote Share", "Committee Chair?", 
                       "Sub-Committee Chair?", 
                       "Years in Congress (log)", 
                       "State Legislature?", "State Legislature Professionalism", 
                       "Average Bills Sponsored (log)", 
                       "Legislative Effectiveness Scores",
                       "Lost Election",
                       "Topic Specialization",
                       "On Power Committee?", "Average Funds Raised",
                       "Cosponsor Centrality")

stargazer(as.data.frame(cor_mat), 
          type = "text", # comment out for tex table
          summary = FALSE,
          covariate.labels = c("Ideological Partisanship", 
                               "Majority", 
                               "Female", "African-American", "Latinx", 
                               "Vote Share", "Chair?", 
                               "Sub-Chair?", 
                               "Years", 
                               "State Legislature?", "Professionalism", 
                               "N Bills Sponsored", 
                               "LES",
                               "Lost",
                               "Specialization",
                               "Power Committee?", "Funds",
                               "Centrality"),
          # uncomment to export table
          #out = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/cor_mat.tex",
          no.space = TRUE)

#---------------
# Figure B2

pol <- filter(pol, cong>104)

tr <- pol %>%
  group_by(cong) %>%
  summarise(rate = n()/535,
            board_rate = sum(seat_next_5yr,na.rm=T)/n())


p <- ggplot(tr, aes(x = cong, y = rate)) +
  geom_line() +
  theme_classic() +
  scale_y_continuous(limits = c(0,.4))  +
  scale_x_continuous( breaks =105:113) +
  labs(x = "Congress", y = "Turnover Rate")

p

# 
# ggsave(p,
#        filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/turnover.eps",
#        device = cairo_ps,
#        width = 8, height = 5)



#---------------
# Table: Board seats of Members of Congress by year

pol$extremist <- ifelse(pol$abs_nom > quantile(pol$abs_nom, 0.74, na.rm = T), 1, 0)

no_miss <- pol %>% 
  drop_na(seat_next_5yr,abs_nom ,
            les , n_bills ,
            chair , sub_chair , av_hhi ,
            seniority ,
            power2 ,
            av_connect ,
            av_fund ,
            state_leg , state_leg_prof , 
            
            lost , chamber ,  party_txt ,  
            majority , female , african_american ,
            latinx , 
            vote_share ,
            state_po,cong)

house_dat <- no_miss %>% filter(chamber == "house")
sen_dat <- no_miss %>% filter(chamber == "senate")

# make tables in house for non-extremists
table(house_dat$last_year[house_dat$extremist == 0], 
      house_dat$seat_next_5yr[house_dat$extremist == 0])

# make tables in house for extremists
table(house_dat$last_year[house_dat$extremist == 1], 
      house_dat$seat_next_5yr[house_dat$extremist == 1])

# make tables in senate for non-extremists
table(sen_dat$last_year[sen_dat$extremist == 0], 
      sen_dat$seat_next_5yr[sen_dat$extremist == 0])

# make tables in senate for extremists
table(sen_dat$last_year[sen_dat$extremist == 1], 
      sen_dat$seat_next_5yr[sen_dat$extremist == 1])


#---------------
# Table: Board seats of Members of Congress in the pre-HLOGA vs. post-HLOGA periods 

table(house_dat$hloga[house_dat$extremist == 0], 
      house_dat$seat_next_5yr[house_dat$extremist == 0])

table(house_dat$hloga[house_dat$extremist == 1], 
      house_dat$seat_next_5yr[house_dat$extremist == 1])


table(sen_dat$hloga[sen_dat$extremist == 0], 
      sen_dat$seat_next_5yr[sen_dat$extremist == 0])

table(sen_dat$hloga[sen_dat$extremist == 1], 
      sen_dat$seat_next_5yr[sen_dat$extremist == 1])


#---------------------
# Full regression table

stargazer(mod, mod,mod_5,mod_5, 
          coef = list(as.vector(coef(mod)),
                      as.vector(c(0,mar_mod_ordered$AME)),
                      as.vector(coef(mod_5)),
                      as.vector(c(0,mar_mod_5_ordered$AME))),
          se = list(as.vector(sqrt(diag(vcov(mod)))),
                    as.vector(c(0,mar_mod_ordered$SE)),
                    as.vector(sqrt(diag(vcov(mod_5)))),
                    as.vector(c(0,mar_mod_5_ordered$SE))),
          omit = c("Constant", "state_po", "factor", "female", 
                   "african_american", "latinx", "party","chamber"),
          covariate.labels = c("Ideological Partisanship",
                               "Legislative Effectiveness Score",
                               "ln Number of Bills", 
                               "Committee Chair", "Subcommittee Chair",
                               "Bill Topic Specialization", "ln Years in Congress",
                               "On Power Committee", "Cosponsor Centrality",
                               "ln Average Campaign Funds",
                               "State Legislature", "State Legislature Professionalism",
                               "Lost Re-Election", "Proportion in Majority", "Average Vote Share"),
          dep.var.labels = c("Immediate Seat", "Seat Within 5 Years"),
          column.labels = c("Log-Odds", "AME", "Log-Odds", "AME"),
          title = "Politician Characteristics and Board Service",
          label = "tab:res",
          type="text", # comment out to get latex table
          #out = "MainTable.tex", # uncomment to export table
          add.lines = list(c("State Fixed Effects?", "Yes", "Yes", "Yes", "Yes"),
                           c("Controls?", "Yes", "Yes", "Yes", "Yes")))


#----------------------
# Standardized regression table


mod_dat<- select(pol, seats_t_ind, seat_next_5yr, 
                 abs_nom , majority , female , african_american ,
                 latinx , 
                 vote_share , chair , sub_chair , seniority ,
                 state_leg , state_leg_prof , 
                 n_bills , 
                 les , lost , chamber , av_hhi , party_txt , power2 , 
                 av_fund, av_connect, state_po)

for(i in c(3:16,18,20:22)){
  
  mod_dat[, i] <- scale(mod_dat[,i])
  
}

sd_mod_5 <- glm(seat_next_5yr ~ 
                  abs_nom + # extremity
                  # expertise
                  les + n_bills +
                  chair + sub_chair + av_hhi +
                  #visibility + connex
                  seniority +
                  power2 +
                  av_connect +
                  av_fund +
                  #other factors and controls
                  state_leg + state_leg_prof + 
                  
                  + lost + chamber +  party_txt +  
                  majority + female + african_american +
                  latinx + 
                  vote_share +
                  state_po , 
                data = mod_dat, family = "binomial")

sd_mod_1 <- glm(seats_t_ind ~ 
                  abs_nom + # extremity
                  # expertise
                  les + n_bills +
                  chair + sub_chair + av_hhi +
                  #visibility + connex
                  seniority +
                  power2 +
                  av_connect +
                  av_fund +
                  #other factors and controls
                  state_leg + state_leg_prof + 
                  
                  + lost + chamber +  party_txt +  
                  majority + female + african_american +
                  latinx + 
                  vote_share +
                  state_po , 
                data = mod_dat, family = "binomial")


stargazer(sd_mod_1, sd_mod_5, 
          omit = c("Constant", "state_po", "factor", "female", 
                   "african_american", "latinx", "party","chamber"),
          covariate.labels = c("Ideological Partisanship",
                               "Legislative Effectiveness Score",
                               "ln Number of Bills", 
                               "Committee Chair", "Subcommittee Chair",
                               "Bill Topic Specialization", "ln Years in Congress",
                               "On Power Committee", "Cosponsor Centrality",
                               "ln Average Campaign Funds",
                               "State Legislature", "State Legislature Professionalism",
                               "Lost Re-Election", "Proportion in Majority", "Average Vote Share"),
          dep.var.labels = c("Immediate Seat", "Seat Within 5 Years"),
          title = "Politician Characteristics and Board Service",
          label = "tab:res",
          type="text", # comment out for tex table
          #out = "StdRegression.tex", # uncomment to export
          add.lines = list(c("State Fixed Effects?", "Yes", "Yes", "Yes", "Yes"),
                           c("Race, Gender, Party, Chamber Controls?", "Yes", "Yes", "Yes", "Yes")))


#----------------------------
# Figure with predicted probabilities (panel A) and probabilities within categories (panel B)

mod <- glm(seats_t_ind ~ abs_nom + majority + female + african_american +
             latinx + 
             vote_share + chair + sub_chair + seniority +
             state_leg + state_leg_prof + 
             n_bills + 
             les + lost + chamber + av_hhi + party_txt + power2 + 
             av_fund + state_po, 
           #log(av_donation+1), 
           data = pol, family = "binomial")

nom_plot <- cplot(mod, "abs_nom", draw=F)

p1<-ggplot(nom_plot, aes(x = xvals, y = yvals)) +
  geom_line() +
  theme_classic() +
  labs(x = "Ideological Partisanship", 
       y = "Predicted Probability of\nBoard Service",
       title = "A: Predicted Probabilities (from logit)")


prop_board <- pol %>%
  mutate(cat_nom = case_when(abs_nom < quantile(abs_nom, probs = .25, na.rm = T) ~ 0,
                             abs_nom > quantile(abs_nom, probs = .75, na.rm = T) ~ 1,
                             TRUE ~ -1)) %>%
  group_by(cat_nom) %>%
  summarise(prob = mean(seats_t_ind, na.rm=T))

p2<-ggplot(prop_board, aes(x = reorder(cat_nom, -prob), y = prob)) +
  #geom_line(group = 1) +
  geom_bar(stat = "identity", fill = "grey", color = "grey") +
  geom_point() +
  scale_x_discrete(labels = c("Moderates", "Loyalists", "Extremists")) +
  labs(x = NULL, y = "Probability of Obtaining an\nImmediate Board Seat",
       title = "B: Probabilities within Categories") +
  theme_classic()


p <- plot_grid(p1, p2)

# print graph
p

# uncomment to export
# ggsave(p,
#        filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/ExtremeViz.eps",
#        width = 8, height = 5)




#------------------------------------------------------------
# Table: "effect of partisanship may be driven by credibility"
#losers vs. retirees

#losers
lost_mod <- glm(seats_t_ind ~ 
                  abs_nom , 
                data = filter(pol, lost==1), family = "binomial")


#retirees 
win_mod <- glm(seats_t_ind ~ 
                 abs_nom , 
               data = filter(pol, lost == 0), family = "binomial")

# interaction with losing
inter_win <- glm(seats_t_ind ~ 
                   abs_nom*lost ,
                 data = pol, family = "binomial")

stargazer(lost_mod, win_mod, inter_win,
                      keep = "abs_nom",
                      type = "text",
                      dep.var.labels = "Immediate Seat",
                      covariate.labels = c("Ideological Extremism",
                                           "Ideological Extremism X Election Loser"),
                      column.labels = c("Election Losers", "Retirees", "Interaction"),
                      add.lines = list(c("Covariates?", "No", "No", "No"),
                                       c("State FE?", "No", "No", "No")),
                      title = c("Effect of Extremism May Be Driven by Credibility"),
                      label = "tab:cred")

#-------------
# with vote share instead of losing indicator

# below 1st quartile
bel_mod <- glm(seats_t_ind ~ 
                 abs_nom , 
               data = filter(pol, vote_share < quantile(pol$vote_share, prob = 0.25, na.rm=T)), family = "binomial")

# above third
ab_mod <- glm(seats_t_ind ~ 
                abs_nom , 
              data = filter(pol, vote_share > quantile(pol$vote_share, prob = 0.75, na.rm=T)), family = "binomial")

#interaction with above 3rd indicator
inter_share <- glm(seats_t_ind ~ 
                     abs_nom*vote_share , 
                   data = pol, family = "binomial")

 stargazer(bel_mod, ab_mod, inter_share,
                        keep = "abs_nom",
                        type = "text",
                        dep.var.labels = "Immediate Seat",
                        covariate.labels = c("Ideological Extremism",
                                             "Ideological Extremism X Vote Share"),
                        column.labels = c("Low Share", "High Share", "Interaction"),
                        add.lines = list(c("Covariates?", "No", "No", "No"),
                                         c("State FE?", "No", "No", "No")),
                        title = c("Effect of Extremism May Be Driven by Credibility"),
                        label = "tab:cred")


#---------------------------------------
# Figure: Distribution of Firm Political Preferences
 
 pol_firm_av <- read_csv(here(this_dir,"data","cfscore_firm_mc.csv"))
 
 # compare distributions
 p<-ggplot(pol_firm_av) +
   geom_histogram(aes(x = cfscore), fill = "red", alpha = .2) +
   geom_histogram(aes(x= recipient.cfscore.dyn), fill = "blue", alpha = .2) +
   theme_classic() +
   annotate(geom = "text", x = -0.5, y = 6.5,
            label = "Politicians")+
   annotate(geom = "text", x = 0.25, y = 8.5,
            label = "Firms") +
   labs(x = "Ideology\n(Campaign Finance Score)", 
        y = "Count")
 
 p
 
 # ggsave("C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/cfscore_hist.eps",
 #        device = cairo_ps,
 #        width = 6, height = 5)
 
#-------------------------------------
# Table: Correlation between Legislator and Firm Ideologies
 
 
 reps <- pol_firm_av %>% filter(party_num == 200)
 
 dems <- pol_firm_av %>% filter(party_num == 100)
 
 dem_mod <- lm(recipient.cfscore.dyn ~ cfscore, data = dems)
 
 rep_mod <- lm(recipient.cfscore.dyn ~ cfscore, data = reps)
 
 
 stargazer(dem_mod, rep_mod, keep.stat = "n",
           type = "text",
           dep.var.labels = "Legislator CFScore",
           covariate.labels = "Firm CFScore",
           title = "Correlation between Legislator and Firm Ideology",
           label = "tab:cfscore",
           column.labels = c("Democrats", "Republicans"))
 #out = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/cfscore_corr.tex")
 
 
 
 
#-------------------------------------
# Figure: Relation Between Legislator and Firm Ideologies.
 
 
 p_dem <- ggplot(dems, aes(x= cfscore, y = recipient.cfscore.dyn)) +
   geom_point() +
   theme_classic() +
   geom_abline(intercept = 0, slope = 1)+
   labs(x = "Firm CFScore", y = "Legislator CFScore",
        title = "A: Democrats")
 
 p_rep <- ggplot(reps, aes(x= cfscore, y = recipient.cfscore.dyn)) +
   geom_point() +
   theme_classic() +
   geom_abline(intercept = 0, slope = 1) +
   labs(x = "Firm CFScore", y = "Legislator CFScore",
        title = "B: Republicans")
 
 
 p <- plot_grid(p_dem, p_rep)
 
 p
 
 # ggsave("C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/cfscore_scatter.eps",
 #        plot =p, 
 #        device = cairo_ps,
 #        width = 10, height = 5)
 
 
 #-----------------------------------------------
 # Results by party
 
 
 mod_inter <- glm(seats_t_ind ~ 
                    abs_nom*dem + # extremity
                    # expertise
                    les + n_bills +
                    chair + sub_chair + av_hhi +
                    #visibility + connex
                    seniority +
                    power2 +
                    av_connect +
                    av_fund +
                    #other factors and controls
                    state_leg + state_leg_prof + 
                    + lost + party_txt +  
                    majority + female + african_american +
                    latinx + 
                    vote_share, 
                  data = pol, family = "binomial")
 

 mod_dem <- glm(seats_t_ind ~ 
                  abs_nom + # extremity
                  # expertise
                  les + n_bills +
                  chair + sub_chair + av_hhi +
                  #visibility + connex
                  seniority +
                  power2 +
                  av_connect +
                  av_fund +
                  #other factors and controls
                  state_leg + state_leg_prof + 
                  + lost + party_txt +  
                  majority + female + african_american +
                  latinx + 
                  vote_share, 
                data = filter(pol, dem == 1), family = "binomial")

 mod_not_dem <- glm(seats_t_ind ~ 
                      abs_nom + # extremity
                      # expertise
                      les + n_bills +
                      chair + sub_chair + av_hhi +
                      #visibility + connex
                      seniority +
                      power2 +
                      av_connect +
                      av_fund +
                      #other factors and controls
                      state_leg + state_leg_prof + 
                      + lost + party_txt +  
                      majority + female + african_american +
                      latinx + 
                      vote_share, 
                    data = filter(pol, dem == 0), family = "binomial")
 
 stargazer(mod_inter, mod_dem, mod_not_dem,
           keep = c("abs_nom"),
           type = "text",
           covariate.labels = c("Ideological Extremism",
                                "Ideological Extremism X Democrat"),
           add.lines = list(c("Controls?", "Yes", "Yes"),
                            c("State FE?", "No", "No")),
           title = "Varying Effects Between Parties",
           #out = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/ByParty.tex",
           label = "tab:dem",
           column.labels = c("Interaction", "Democrats", "R+I"),
           dep.var.labels = "Immediate Seat")
 
 
 
#-----------------------------------------------------------
# Effects by chamber

mod_inter <- glm(seats_t_ind ~ 
                   abs_nom*chamber + # extremity
                   # expertise
                   les + n_bills +
                   chair + sub_chair + av_hhi +
                   #visibility + connex
                   seniority +
                   power2 +
                   av_connect +
                   av_fund +
                   #other factors and controls
                   state_leg + state_leg_prof + 
                   + lost + party_txt +  
                   majority + female + african_american +
                   latinx + 
                   vote_share, 
                 data = pol, family = "binomial")



mod_house <- glm(seats_t_ind ~ 
                   abs_nom + # extremity
                   # expertise
                   les + n_bills +
                   chair + sub_chair + av_hhi +
                   #visibility + connex
                   seniority +
                   power2 +
                   av_connect +
                   av_fund +
                   #other factors and controls
                   state_leg + state_leg_prof + 
                   + lost + party_txt +  
                   majority + female + african_american +
                   latinx + 
                   vote_share, 
                 data = filter(pol, chamber == "house"), family = "binomial")


mod_sen <- glm(seats_t_ind ~ 
                 abs_nom + # extremity
                 # expertise
                 les + n_bills +
                 chair + sub_chair + av_hhi +
                 #visibility + connex
                 seniority +
                 power2 +
                 av_connect +
                 av_fund +
                 #other factors and controls
                 state_leg + state_leg_prof + 
                 + lost + party_txt +  
                 majority + female + african_american +
                 latinx + 
                 vote_share, 
               data = filter(pol, chamber == "senate"), family = "binomial")

stargazer(mod_inter, mod_house, mod_sen,
          keep = c("abs_nom"),
          type = "text",
          covariate.labels = c("Ideological Extremism",
                               "Ideological Extremism X Senate"),
          add.lines = list(c("Controls?", "Yes", "Yes"),
                           c("State FE?", "No", "No")),
          title = "Varying Effects Between Chambers",
          #out = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/ByChamber.tex",
          label = "tab:chamber",
          column.labels = c("Interaction", "House", "Senate"),
          dep.var.labels = "Immediate Seat")


#------------------------------------
# Effects over time


dv <- data.frame(seats_t_ind = pol$seats_t_ind, 
                 seats_t1_ind = pol$seats_t1_ind, 
                 seats_t3_ind = pol$seats_t3_ind, 
                 seats_t5_ind = pol$seats_t5_ind)

time_res <- list()

for(i in 1:ncol(dv)){
  
  crnt_mod <- glm(dv[,i] ~ abs_nom + # extremity
                    # expertise
                    les + n_bills +
                    chair + sub_chair + av_hhi +
                    #visibility + connex
                    seniority +
                    power2 +
                    av_connect +
                    av_fund +
                    #other factors and controls
                    state_leg + state_leg_prof + 
                    + lost + chamber +  party_txt +  
                    majority + female + african_american +
                    latinx + 
                    vote_share +
                    state_po, 
                  data = pol, family = "binomial")
  crnt_mar <- summary(margins(crnt_mod, variables = "abs_nom")) 
  
  
  #crnt_mod <- crnt_mar %>% filter(term == "abs_nom")
  crnt_mar$dv <- names(dv)[i]
  
  time_res[[i]] <- crnt_mar 
}

var_names <- labeller(term = 
                        c("abs_nom" = "Ideological Partisanship"))

time_res <- do.call("rbind.data.frame", time_res)

time_res$time <- 1:4

p_time <- ggplot(time_res, aes(x = time, y = AME)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower,
                    ymax = upper, width = 0)) +
  geom_line() +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Years Since Leaving Office",
       y= "Average Marginal Effect of\nIdeological Partisanship on Board Service") +
  scale_x_continuous(breaks = 1:4,
                     labels = c("0", "1", "3", "5"))




p_time

# uncomment to export
#ggsave(filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/time_trace.eps",
#       width = 6.5, height = 5)


#----------------------------------
# table with event history model
# i.e. year FE model for robustness
year_fe_mod <- glm(seats_t_ind ~ 
                     abs_nom + factor(cong), 
                   data = pol, family = "binomial")

stargazer(year_fe_mod,
          covariate.labels = "Ideological Extremism",
          title = "Beck et al Event History Estimates",
          keep.stat = "n", keep = "abs_nom", 
          type = "text",
          dep.var.labels = "Immediate Board Seat")
          #out = "EventHistory.tex")



#----------------------------------
# Sensitivity to excluding observations

# note these models cannot contain state fixed effects
excl_df <- pol %>% drop_na(seats_t_ind, 
                       abs_nom,les , n_bills ,
                       chair , sub_chair , av_hhi ,
                       #visibility , connex
                       seniority ,
                       power2 ,
                       av_connect ,
                       av_fund ,
                       #other factors and controls
                       state_leg , state_leg_prof , 
                       lost , chamber ,  party_txt ,  
                       majority , female , african_american ,
                       latinx , 
                       vote_share)

reg_form <- formula(seats_t_ind ~ 
                      abs_nom + # extremity
                      # expertise
                      les + n_bills +
                      chair + sub_chair + av_hhi +
                      #visibility + connex
                      seniority +
                      power2 +
                      av_connect +
                      av_fund +
                      #other factors and controls
                      state_leg + state_leg_prof + 
                      
                      + lost + chamber +  party_txt +  
                      majority + female + african_american +
                      latinx + 
                      vote_share)
#state_po)
biv_mod <- lm(formula=reg_form, x=TRUE, y=TRUE, 
              data = excl_df)
#summary(biv_mod)

model_grads <- ComputeModelInfluence(biv_mod) %>%
  AppendTargetRegressorInfluence("abs_nom")


signals <- GetInferenceSignals(model_grads)

# Predict the changes, re-run the model at the left-out points, and
# compare the two.

preds <- PredictForSignals(signals, model_grads)
reruns <- RerunForSignals(signals, model_grads)
reruns_df <- GetSignalsAndRerunsDataframe(signals, reruns, model_grads)
preds_df <- GetSignalsAndRerunsDataframe(signals, preds, model_grads)

summary_df <-
  rbind(reruns_df %>% mutate(method="rerun"),
        preds_df %>% mutate(method="prediction")) %>%
  pivot_wider(names_from=method, values_from=value)

summary_df <- summary_df[c(1,5,9),]

summary_df <- summary_df %>% select(description, prop_drop)
summary_df$prop_drop <- round(summary_df$prop_drop, digits = 3)


stargazer(summary_df, summary = FALSE,
          covariate.labels = c("", "Target Change", "Proportion Dropped"),
          title = "Broderick et al (2020) Sensitivity Analysis",
          label = "sens_out",
          type = "text")
          #out = "sensitivity_to_exclusion.tex")

signal <- signals[["abs_nom"]][["both"]]
excl_df$drop <- GetWeightVector(signal$apip$inds, nrow(excl_df), bool=TRUE, invert=TRUE)
excl_df$infl <- signal$qoi$infl


sens_plot <- grid.arrange(
  ggplot(excl_df) +
    geom_point(aes(x=abs_nom, y=infl, color=drop)) +
    labs(x = "Ideological Partisanship", y = "Influence") +
    theme_bw(),
  ggplot(excl_df) +
    geom_point(aes(x=abs_nom, y=seats_t_ind, color=drop))+
    labs(x = "Ideological Partisanship", y = "Immediate Board Seat")+
    theme_bw(),
  ncol=2
)

#uncomment to export graph
# ggsave(plot = sens_plot,
#        filename = "SensitivityOutliers.eps",
#        width = 10, heigh = 5)
# 


#--------------------------------------------
# Robustness to the choice of controls

mod <- glm(seats_t_ind ~ 
             abs_nom , 
           data = pol, family = "binomial")

glm_logit <- function(formula, data) {
  glm(formula = formula, 
      data = data, 
      family = "binomial")
}

results <- run_specs(df = pol, 
                     y = "seats_t_ind", 
                     x = "abs_nom", 
                     model = c("glm_logit"), 
                     controls = c("les", "n_bills",
                                  "chair", "sub_chair", "av_hhi",
                                  "seniority",
                                  "power2",
                                  "av_connect",
                                  "av_fund",
                                  "state_leg", "state_leg_prof", 
                                  "lost", "chamber",  "party_txt",  
                                  "majority", "female", "african_american",
                                  "latinx", 
                                  "vote_share",
                                  "state_po"))

progressive_plot <- plot_specs(results)

progressive_plot

# ggsave(plot = progressive_plot,
#        filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/progressive_controls.eps",
#        device = cairo_ps,
#        width = 8,
#        height = 10)



#------------------------------------------------------------------
# Robustness of the difference-in-differences results to grouping


pol <- as.data.frame(pol)

alt_bw <- c(0.977,
            0.5*0.977,
            0.4*0.977,
            0.3*0.977,
            0.2*0.977,
            0.1*0.977)


kern_func <- function(this_bw){
  
  
  eks_kern <- interflex(estimator = "kernel",
                        data = pol, 
                        Y = "seats_t_ind",
                        D = "did", X = "abs_nom",
                        Z = c("sen", "hloga"),
                        CI=FALSE, bw= this_bw,
                        parallel = T, cores = 8,nboots=200,
                        cl="state_po",
                        na.rm=T)

  
  res <- as.data.frame(eks_kern$est.kernel)
  res$bw <- this_bw
  return(res)
  
}


four_bw <- map_df(alt_bw,
                  ~ kern_func(.))
names(four_bw) <- c("abs_nom", "ME")

four_bw$bw_txt <- rep(c("A: Optimal BW (0.98)", "B: BW = 0.49",
                        "C: BW = 0.38", "D: BW = 0.29",
                        "E: BW = 0.19", "F: BW = 0.1"), each = 50)

p<- ggplot(four_bw, aes(x = abs_nom, y = ME)) +
  geom_line() +
  facet_wrap(~bw_txt) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Ideological Partisanship\n(Absolute Distance to Chamber-Period Median of DW-NOMINATE)",
       y = "ATT of HLOGA on Board Service")

p

# ggsave(plot = p,
#        filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/different_bw.eps",
#        device = cairo_ps,
#        width = 10,
#        height = 6)


#--------------------------------------------------
# Table with HLOGA DiD results


# extract estimates for graph
df_tab <- filter(df_ex, spec == "HLOGA") %>% select(-c("adj_x0", "sd", "spec"))

names(df_ex) <- c("Extremism", "Diff-in-Diff", "Lower CI", "Upper CI")
rownames(df_tab) <- NULL

stargazer(df_tab, 
          type = "text",
          summary = FALSE,
          title = "HLOGA Diff-in-Diff Estimates",
          label = "tab:DD_table")
          #out = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/HLOGA_DD_table.tex")


#--------------------------------------------------------
# Table: Strength of Partisanship and Think Tank and Media Jobs


ThinkMed_df <- pol %>% drop_na(seats_t_ind , 
          abs_nom , # extremity
          # expertise
          les , n_bills ,
          chair , sub_chair , av_hhi ,
          #visibility , connex
          seniority ,
          power2 ,
          av_connect ,
          av_fund ,
          #other factors and controls
          state_leg , state_leg_prof , 
          lost , chamber ,  party_txt ,  
          majority , female , african_american ,
          latinx , 
          vote_share,
          state_po , cong)

think <- read_csv(here(this_dir,"data","ThinkTankData.csv"))

think$think_tank <- ifelse(is.na(think$think_tank)==T, 0, think$think_tank)

think$five_year <- think$start_year - think$last_year

think <- think %>%
  filter(think_tank == 1 & five_year >-1 & five_year < 5) %>%
  group_by(id_unique) %>%
  summarise(think_tank = mean(think_tank,na.rm=T)) %>%
  mutate(think_tank = ifelse(think_tank>0,1,0))



media <- read_excel(here(this_dir, "data","MediaData.xlsx"))


media$five_year <- media$start_year_nodate - media$last_year

media <- media %>%
  filter(Talking_head == 1 & five_year >-1 & five_year < 5) %>%
  group_by(id_unique) %>%
  summarise(Talking_head = mean(Talking_head,na.rm=T)) %>%
  mutate(Talking_head = ifelse(Talking_head>0,1,0))


# merge in
ThinkMed_df <- left_join(ThinkMed_df, think, by = c("id" = "id_unique"))
ThinkMed_df <- left_join(ThinkMed_df, media, by = c("id" = "id_unique"))


ThinkMed_df$think_tank <- ifelse(is.na(ThinkMed_df$think_tank) == T, 0, ThinkMed_df$think_tank)
ThinkMed_df$Talking_head <- ifelse(is.na(ThinkMed_df$Talking_head ) == T, 0, ThinkMed_df$Talking_head )

#----------------------
# results

think_mod <- glm(think_tank ~ abs_nom,
                 data = ThinkMed_df, family = "binomial")
summary(think_mod)
think_mar <- summary(margins(think_mod))


media_mod <- glm(Talking_head ~ abs_nom,
                 data = ThinkMed_df, family = "binomial")
summary(media_mod)
media_mar <- summary(margins(media_mod))


stargazer(think_mod,media_mod,
          keep.stat = "n",
          type = "text",
          add.lines = list(c("AME",  round(think_mar$AME,digits = 3),
                             round(media_mar$AME,digits = 3))),
          covariate.labels = "Ideological Extremism", 
          dep.var.labels = c("Think Tank Job", "Media Job"),
          title = "Extremism and Think Tank and Media Jobs",
          label = "tab:jobs")
          #out = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/logit_jobs.tex")

#---------------------------------------------------------------------------------------
# Figure: The HLOGA Diff-in-Diff for Jobs in Media an Think Tanks.

ThinkMed_df <- as.data.frame(ThinkMed_df)
did_think <- interflex(estimator = "binning",
                       data = ThinkMed_df, Y = "think_tank",
                       D = "did", X = "abs_nom",cutoffs = quantile(ThinkMed_df$av_connect, probs = c(0.25,0.5,0.75),na.rm=T),
                       Z = c("sen", "hloga"), vartype="delta",
                       cl="state_po",vcov.type = "cluster",
                       FE = "cong", na.rm=T, full.moderate = F)


att_think <-ggplot(as.data.frame(did_think$est.bin$`1`),
                   aes(x = x0, y = coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = coef - 1.96*sd,
                    ymax = coef + 1.96*sd), width = 0, colour = "grey") +
  geom_errorbar(aes(ymin = coef - 1.645*sd,
                    ymax = coef + 1.645*sd), width = 0) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Ideological Partisanship\n(Absolute Distance to Chamber-Period Median of DW-NOMINATE)",
       y = "ATT of HLOGA on Board Service",
       title = "A: Think Tank Jobs")


did_media <- interflex(estimator = "binning",
                       data = ThinkMed_df, Y = "Talking_head",
                       D = "did", X = "abs_nom",cutoffs = quantile(ThinkMed_df$av_connect, probs = c(0.25,0.5,0.75),na.rm=T),
                       Z = c("sen", "hloga"), vartype="delta",
                       cl="state_po",vcov.type = "cluster",
                       FE = "cong", na.rm=T, full.moderate = F)


att_media <-ggplot(as.data.frame(did_media$est.bin$`1`),
                   aes(x = x0, y = coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = coef - 1.96*sd,
                    ymax = coef + 1.96*sd), width = 0, colour = "grey") +
  geom_errorbar(aes(ymin = coef - 1.645*sd,
                    ymax = coef + 1.645*sd), width = 0) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Ideological Partisanship\n(Absolute Distance to Chamber-Period Median of DW-NOMINATE)",
       y = "ATT of HLOGA on Board Service",
       title = "B: Talking Head Jobs")


p_jobs <- plot_grid(att_think, att_media)
p_jobs

# ggsave(plot = p_jobs,
#        filename = "C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/alternative_jobs.eps",
#        device = cairo_ps,
#        width = 9, height = 5)


#---------------------------------------------------------------------
# Figure: The HLOGA Decreases Voluntary Retirements among Extremist Senators.


mc_nom2 <- read_csv(here(this_dir,"data","voluntary_retirement_repli.csv"))


set.seed(1010221102)
mod <- felm(leave ~ extremist*hloga*chamber
            |0|0|id, data = mc_nom2,
            Nboot = 500,
            nostat=structure(FALSE, boot= T),
            bootcluster = factor(mc_nom2$id))

summary(mod)

reps <- mod$boot
reps <- as.data.frame(t(reps))

#-----
set.seed(101020)

sen_dat <- filter(mc_nom2, chamber == "senate")
sen_mod <- felm(leave ~ extremist*hloga
                |0|0|id, data = sen_dat,
                Nboot = 500,
                nostat=structure(FALSE, boot= T),
                bootcluster = factor(sen_dat$id))
summary(sen_mod)

sen_reps <- sen_mod$boot
sen_reps <- as.data.frame(t(sen_reps))
quantile(sen_reps$`extremist:hloga`, probs =c(0.025,0.05,0.95,0.975), na.rm=T)

house_dat <- filter(mc_nom2, chamber == "house")

set.seed(101020)

house_mod <- felm(leave ~ extremist*hloga
                  |0|0|id, data = house_dat,            
                  Nboot = 500,
                  nostat=structure(FALSE, boot= T),
                  bootcluster = factor(house_dat$id))
summary(house_mod)

house_reps <- house_mod$boot
house_reps <- as.data.frame(t(house_reps))


triple <- data.frame(est = c(coef(mod)[8],
                             coef(sen_mod)[4],
                             coef(house_mod)[4]),
                     rbind(quantile(reps$`extremist:hloga:chambersenate`, probs =c(0.025,0.05,0.95,0.975), na.rm=T),
                           quantile(sen_reps$`extremist:hloga`, probs =c(0.025,0.05,0.95,0.975), na.rm=T),
                           quantile(house_reps$`extremist:hloga`, probs =c(0.025,0.05,0.95,0.975))),
                     spec = c("Triple Difference",
                              "Senate DiD",
                              "House DiD"))

names(triple)[2:5] <- c("lwr95", "lwr90", "upr90", "upr95")

p<-ggplot(triple, aes(x = est, y = spec)) +
  geom_point() +
  geom_errorbarh(aes(xmin = lwr95, xmax = upr95), colour = "grey",
                 height = 0) +
  geom_errorbarh(aes(xmin = lwr90, xmax = upr90), colour = "black",
                 height=0) +
  theme_classic() +
  geom_vline(xintercept = 0, lty = 3) +
  labs(x = "Estimated ATT", y = NULL)

p

# ggsave("C:/Users/bcke.egb/Dropbox (CBS)/Politician-Directors/latex/images/voluntary_retirement.eps",
#        plot =p, 
#        device = cairo_ps,
#        width = 6, height = 5)


#----------------------------------------------------------------------------
# Table: Has HLOGA Affected Retirement Overall?

mod_tot <- felm(leave~ hloga*chamber
                |0|0|id, data = mc_nom2,
                Nboot = 500,
                nostat=structure(FALSE, boot= T),
                bootcluster = factor(mc_nom2$id))


stargazer(mod_tot, keep.stat = "n",
          type = "text",
          covariate.labels = c("Post-HLOGA",
                               "Senator", 
                               "Post-HLOGA X Senator"),
          dep.var.labels = "Voluntary Retirement",
          title = "Has HLOGA Affected Retirement Overall?",
          label = "tab:HLOGA_trend")

sink()








